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With quantum computers being out of reach for now, quantum simulators are the alternative 
devices for efficient and more exact simulation of problems that are challenging on conventional 
computers. Quantum simulators are classified into analog and digital, with the possibility of con¬ 
structing “hybrid” simulators by combining both techniques. In this paper, we focus on analog 
quantum simulators of open quantum systems and address the limit that they can beat classical 
computers. In particular, as an example, we discuss simulation of the chlorosome light-harvesting 
antenna from green sulfur bacteria with over 250 phonon modes coupled to each electronic state. 
Furthermore, we propose physical setups that can be used to reproduce the quantum dynamics of a 
standard and multiple-mode Holstein model. The proposed scheme is based on currently available 
technology of superconducting circuits consist of flux qubits and quantum oscillators. 

PACS numbers: 71.38.Ht; 03.67.Ac; 74.25.Ld; 74.81.Fa. 


There is a growing interest in understanding the dy¬ 
namics of open quantum systems, particularly, when a 
particle is coupled to a vibrational environment. Such sit¬ 
uations arise in quantum chemistry and condensed mat¬ 
ter physics, for example, in photosynthetic complexes or 
molecular aggregates. Thus, a detailed study of the dy¬ 
namics of electron-phonon interaction becomes desirable. 
Although many analytical and numerical methods have 
been applied to this problem [1-11], their applicability is 
often limited by the number of the phonon modes cou¬ 
pled to the electronic states or to a particular investiga¬ 
tion (e.g. low-lying excited polaron) and a specific pa¬ 
rameter regime. The resources required for most of the 
classical computational methods increase exponentially 
with the number of particles in the simulation and it is 
challenging to simulate the dynamics of open quantum 
systems on conventional computers, even using modern 
parallel processing units [12-14]. The situation becomes 
even much more challenging for complex open quantum 
systems with structured environments. As yet, only small 
model systems have been studied theoretically with crude 
approximations to the system-bath dynamics, see for ex¬ 
ample [15, 16]. Numerically exact solution can be ob¬ 
tained for only small systems (< 20 sites) with restric¬ 
tions on the bath modes [12, 14, 17-20]. In Figure 1, 
we estimate the upper limit for simulating such complex 
open quantum systems with current computational re¬ 
sources on conventional computers. The horizontal axis 
indicates the system size (number of the particles or sites) 
that can be simulated while the vertical axis indicates 
number of the peaks in the spectral density that could 
be considered in this simulation, see Ref. [14] for more 
computational details. Note that “peaks” here refers to 
Drude-Lorentz peaks in the spectral density [14] which 
should not be confused with the number of phonon modes 
in the Hamiltonian. Each of these peaks in the spectral 
density may include several phonon modes. 
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Figure 1. The grey area shows the estimated treatable 
system sizes for the simulation of Frenkel exciton Hamiltoni¬ 
ans using current classical supercomputing resources. There 
is a trade off between the complexity of the spectral den¬ 
sity and the system size that denotes the classically-feasible 
area. Three photosynthetic systems are shown: The Fenna- 
Mathews Olson (FMO) complex of Green-Sulfur Bacteria, the 
Light Harvesting I and II complexes of Purple Bacteria and 
Photosystem II of higher plants. The simulation has been per¬ 
formed using the hierarchical equations of motion (HEOM) 
approach on 64 AMD Opteron cores employing a total of 250 
GB of RAM. 


In this work, we propose analog quantum devices [21- 
23] to mimic the dynamics of complex open quantum 
systems and demonstrate that they can be constructed 
using present-day technology of superconducting circuits 
and outperform current classical computational meth¬ 
ods. With such quantum simulators, one can perform 
more extensive investigation including exciton transport, 
spectral density, absorption spectra as well as wide range 
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of parameters and thereby a more detailed understand¬ 
ing of the problems. Furthermore, our proposed quan¬ 
tum simulators occupy a wide region in the plot shown 
in Figure 1. Similar ideas for simulating Holstein po- 
larons based on polar molecules trapped in an optical 
lattice [24, 25], Rydberg states of cold atoms and ions 
[26], trapped ions [27, 28] and superconducting circuit 
quantum electrodynamics (QED) [29, 30] have been pur¬ 
sued earlier. However, the main focus of this paper is 
emulating the dynamics of multiple-mode Holstein mod¬ 
els at finite temperature - with application in open quan¬ 
tum systems with complex environments - which has not 
been addressed in any of the above mentioned references. 
Moreover, it is worthwhile to study an alternative set-up, 
since different experimental realizations carry distinct ad¬ 
vantages and drawbacks. 

Standard Holstein model. We first focus on simu¬ 
lating an electron-phonon model which describes the in¬ 
teraction of a single electron on a ID finite lattice with 
one vibrational mode per lattice site, namely the Holstein 
model: 


Had = l + H ph + H e i_ph . (1) 

The first term of the above Hamiltonian (electronic 
term) is given by H e i = En=i V n ( a l a n +1 + 4+i a n) 
with V n being the strength of the nearest-neighbor cou¬ 
plings, a\ (a n ) being the creation (annihilation) op¬ 
erators of the electron and N being the number of 
sites. The phonon Hamiltonian is H p h = Yln=i 
with uj n being the frequency of the phonon mode cou¬ 
pled to the n-th lattice site and b^ n (b n ) being the cre¬ 
ation (annihilation) operators of the phonon. The last 
term in Eq. (1) describes the electron-phonon coupling 

He l—ph = En=i K ri a\a n (&* + b n ) with K n being the cou- 
pling strength between the electron and phonon at lat¬ 
tice site n. Using the Jordan-Wigner transformation, the 
Hamiltonian HhoI can be rewritten in terms of the Pauli 
a operators, 

ffHol = 2 £ v n K<x" +1 + <rX +1 ) + 

n= 1 
N 

+ ^2 [ K n ct" (bi + b n ) + hio n b ] n b n ] . ( 2 ) 

n= 1 

In order to reproduce the quantum dynamics of the 
open system given by the above Hamiltonian, let us 
consider a chain of N gradiometric flux qubits [23, 31] 
with tunable coup lings [32, 33] and a single LC- 

oscillator coupled to each qubit, as shown in Figure 2. 
The Hamiltonian of a single flux qubit in the bare basis, 
the quantum states with magnetic flux pointing up |t) 
and down ||), is given by = (Si a l z + cr*)/2 [34], 

where Si is the energy bias between |t) and ||), A i is 
the tunnel splitting between the two states and i labels 
the position of the qubit in the chain. Note that Si can 
be tuned to zero to neglect the term Si(j l z and therefore 



Figure 2. Superconducting quantum circuit diagram of the 
proposed quantum simulator for the Holstein model. The 
qubit states are encoded in the quantized circulating current 
of the qubit loop. The red crosses denote Josephson junc¬ 
tions. The gradiometric flux qubits are coupled with a tun¬ 
able cr^cr^-coupling. Each of the qubits is independently cou¬ 
pled to a quantum LC-oscillator to simulate the vibrational 
environment. 


be at the optimal operating point [35] of the flux qubit, 
which is the most common case in current experiments. 
The coupling between two nearest-neighbor qubits in 
the bare basis is given by H l couv = ^(A^ +1 ) ct* +1 , 

where A ^ +1 is the (tunable) tunnel splitting of the 
coupler qubit (smaller qubits in Figure 2 , see Ref. [23] 
for more details). The coupling of a quantum LC- 
oscillator to the smaller loop of a flux qubit, as shown 

in Figure 2 , is given by H^_ osc = r]i a x (c\ Hb c^j with 

c\ ( Ci ) being the creation (annihilation) operator of 
the oscillator coupled to the i-th qubit and r]i being 
the coupling strength. Finally, the Hamiltonian of a 
single oscillator is H l osc = c[ci with uj[ being the 
transition frequency of the oscillator. Rewriting the 
above Hamiltonians in the energy eigenbasis of the qubit 
| =b) = (||) zb |t)) /y/2 converts the operators a x —)> o\ 
and a\ <x * +1 ->■ < 7 * <rj , +1 « (<r* a i + 1 + <7* cr* +1 ) / 2 in the 
rotating wave approximation (neglecting strongly 
off-resonant couplings). Then the total Hamil¬ 
tonian of the superconducting circuit proposed 
to emulate the dynamics of the Holstein model 

Hsim = EZl { H \ + #coup + #q-osc + #osc} in the new 

basis is given by 


N-l 


Hsim « \ J2 9i( A ii+ 1) K < +1 + < +1 ] 

-I(: 

i= l ^ 


r a i , , ( 

t \ 

( -y < 7 * + ffc < 7 , (< 

4 +a) 


(3) 


Here we can assume identical flux qubits, similar to the 
most common case considered in the literature of the Hol¬ 
stein model. Then the term A^/2cr* in the above 

Hamiltonian leads to a global phase. Comparison of the 
Hamiltonians ( 2 ) and (3), demonstrates that a chain of 
coupled flux qubits with a single quantum LC-oscillator 
coupled to each qubit can simulate the same dynamics of 
the Holstein model with gi( A^ +1 ), 77 ^, u)[ corresponding 
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to V n , tin, u n , respectively. Interestingly, for supercon¬ 
ducting flux qubits, the couplings ^(AF +1 ) and rji are 
tunable. The implementable range of AF +1 ) is in the 
range of approximately zero to 1 [GHz] [36] . rji can be in 
the <10 [GHz] range depending on the frequency of the 
resonator. The tunability and wide implementable range 
of these parameters makes it possible to study different 
parameter regimes of interest (strong coupling rji 
weak coupling rji <C gi and intermediate rji ~ regimes) 
using the proposed quantum simulator. 

Preparation of the qubits in their ground state is 
straightforward: One needs to allow them to relax as 
close as possible to their ground state by cryostatic cool¬ 
ing. Subsequently, the qubits can be initialized by flux 
control in the appropriate initial states for the simula¬ 
tion. The excitation of a qubit is undemanding to achieve 
with the application of a resonant microwave excitation 
(7r-pulse) carried by a microwave line which is connected 
to the respective qubit. This technique has been used 
extensively, e.g., for the observation of Rabi oscillations 
in a flux qubit [34, 37]. After some evolution time the 
populations of the qubit states are measured. 

Temperature. Note that the standard Holstein 
model discussed above does not contain temperature, 
however, the superconducting circuit (as a real physical 
system) is at finite temperature T s i m . Currently, a su¬ 
perconducting circuit can be refrigerated down to a very 
low temperature, around 10 [mK]^0.2 [GHz], and the 
flux qubits can be even cooled down far below 10 [mK] 
using active microwave cooling [38] . Although the quan¬ 
tum simulator being at finite temperature seems to be a 
disadvantage, we will see in the following that the easy 
tunability over T s i m allows one to investigate the physi¬ 
cally relevant case of a finite temperature Holstein model 
over a wide range of temperatures. To this end, we will 
generalize the Standard Holstein Hamiltonian. 

Generalized Holstein model The Hamiltonian of 
multi-mode Holstein model is given by 


N N 

tfgen =2 EE J » (« + «) 

n= 1 m= 1 
n^m 


N 






where k labels the vibrational modes coupled to the site 
n with frequency u n k, K n k = hUnkVRnk is the coupling 
of the electronic excitation of the site n to the vibra¬ 
tional mode k with R nk being the dimensionless Huang- 
Rhys factor (electron-phonon coupling constant) [39], 
and constant C n = e n + huj nk Rnk + D n with e n being 
the electronic transition energy, and D n being the gas-to- 
crystal shift of the transition energy due to nonresonant 
forces [39, 40]. Since a constant energy offset does not 
alter the dynamics, we ignore C n in our approach. Now 
each site couples to a set of oscillators with frequencies 







Figure 3. Representation of a single flux qubit coupled to 
quantum LC-oscillators. (a) Many single resonators are di¬ 
rectly coupled to the qubit, (b) Using a linear-algebraic bath 
transformation [42], the set of independent resonators (di¬ 
rectly coupled to the qubit) are transformed into a set of 
weakly-coupled multiple parallel chain of resonators. 


uj n k and corresponding couplings K n k- We have also gen¬ 
eralized interactions J nrn between arbitrary sites. The 
dynamics of a multiple-mode Holstein model can be re¬ 
produced by a similar superconducting circuit shown in 
Figure 2 with additional quantum LC-oscillators coupled 
to each flux qubit, see Figure 3 (a). 

The experimental implementation of such a quantum 
simulator can face challenges due to the current con¬ 
straints in the realizable superconducting circuits. The 
number of quantum LC-oscillators, that are directly cou¬ 
pled to a qubit is limited by the physical size of the su¬ 
perconducting qubits. Moreover, the coupling strength of 
the qubit to the quantum oscillator is limited and should 
not exceed a certain percentage of the frequency of the 
oscillator [41]. A simple numerical formula for the cou¬ 
pling K nk is given by, 


^nk 
Unk 



5.48/3 nfc /y / Zf 
50 [nA] \ 100 [Ohms] 


1/2 


Ujik 

27r [GHz] 


.(5) 


Here /3 n k is the inductive division ratio (flux of the k- 
th oscillator coupled to the n-th qubit is f3 n k times of 
the qubit flux). This parameter needs to be far be¬ 
low 1 to avoid hybridizing the qubit with the resonator. 
Z™ k is the oscillator impedance and has to be well be¬ 
low the impedance of free space (not much higher than 
100 [Ohms]), in order to maintain high quality factors 
for the resonators. I™ k is the effective persistent current 
of the DC SQUID loop, which is the linear slope of the 
qubit energy splitting with respect to DC SQUID flux. 
In principle, I™ k can be made large, but this also can 
cause linear flux sensitivity of the qubit energy. 

These challenges can be addressed and resolved by a 
linear-algebraic bath transformation [42] that we have 
proposed recently. Based on a simple linear algebraic 
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approach, the set of independent LC oscillators directly 
coupled to a qubit, Figure 3 (a), can be transformed into 
a set of weakly-coupled multiple parallel chain of oscilla¬ 
tors, see Figure 3 (b). This transformation can dramati¬ 
cally reduce the number of the oscillators that are directly 
coupled to the qubit as well as the coupling strength of 
the quantum oscillators to the qubit. To specify the num¬ 
ber of the required resonators and their parameters and 
to feature outrunning classical algorithms with our pro¬ 
posed approach, as an example, here we study the feasi¬ 
bility and provide an outlook for the emulation of the dy¬ 
namics of the chlorosome light-harvesting antenna from 
green sulfur bacteria. 

Chlorosome light-harvesting antennae. The 

green sulfur bacteria lives in a deep sea where only a few 
hundred photons per second arrive at a bacterium [43]. 
Therefore, they should be able to transfer the photon 
energy efficiently, rapidly and robustly to the reaction 
center to generate the electro-chemical potential energy 
gradient and exploit it in the photosynthetic metabolic 
cycle. Compared with other light-harvesting species, the 
chlorosome has the unique feature that it is composed 
of 200-250 thousands bacteriochlorophyll molecules or¬ 
ganized into supramolecular assembly [44, 45]. How the 
quantum dynamics helps the excitation energy transfer 
within this giant molecular aggregate is an interesting 
question and has attracted many research groups, see 
the references cited in Refs. [15, 16, 44]. The dynam¬ 
ics of the chlorosome can be given by the multi-mode 
Holstein model in Eq.(4). The structure model has been 
proposed from experiments [45] and studied theoretically 
by some of the authors [15, 16] with a crude stochastic 
quantum propagating model and its spectral density is 
demonstrated in Figure 4. 

As discussed above, the dynamics of this system can 
be emulated by a chain of superconducting qubits and 
253 quantum LC-oscillators coupled to each qubit. The 
size of each flux qubit is around few ten of microns and 
there is no enough physical space to couple it directly to 
253 resonators. However, we can reduce number of the 
resonators that are directly coupled to the qubit by using 
the linear-algebraic bath transformation [42] . This trans¬ 
formation mixes resonator modes with different frequen¬ 
cies to distributes 253 modes to, for example, a set of 6 
parallel chains of quantum resonators, Figure 3 (b), with 
each of the chains having the most of 43 coupled oscil¬ 
lators. In addition to reducing number of the resonators 
that are directly coupled to the qubit, this mapping will 
also reduce the coupling strength of the qubit to the pri¬ 
mary oscillator modes (the first oscillators in the chains 
that are directly coupled to the qubit). 

Note that the parameters of the superconduct¬ 
ing simulator are temperature dependent. To 
account for finite temperature, we first trans¬ 
form the spectral density given in Figure 4 using 
C(c<j, T) = {1 + coth [hu/ (2kBT)]} J a (cj), where the 
subscript “A” denotes the antisymmetric spectral den¬ 
sity J a (uj) = J(uj ) if uj > 0; and J a (cj) = — J{—uS) if 


10000 



Frequency [cnr 1 ] 


Figure 4. Spectral density of the electron-phonon coupling 
of bacteriochlorophyll molecules in the chlorosome antenna 
of green-sulfur bacteria with 253 phonon modes. The spec¬ 
tral density of the phonon bath was obtained from a quan¬ 
tum mechanics/molecular mechanics (QM/MM) simulation 
with time-dependent density function theory in Ref. [15]. An 
experiment ally-resolved [45] structure of the chlorosome is 
shown in the inset. 


uj < 0 , see [23] for more details. Since the chlorosome is 
at room temperature, T c h=300 [K], and the supercon¬ 
ducting circuit can be considered at T sim =10 [mK], then 
all the parameters of the quantum simulator need to be 
rescaled accordingly; Aim = (T sim /T ch ) Ah, with Aim 
and A c h indicating any parameter of the quantum simu¬ 
lator and chlorosome, respectively. Then after rescaling, 
we perform the linear-algebraic bath transformation. 
With this procedure, the coupling strengths between the 
qubit and the oscillators that are directly coupled to the 
qubit, Figure 3 (b), need to be around 150 - 210 [MHz]. 
The coupling between the oscillators in the chains are 
around 100 - 560 [MHz], the required frequencies for 
the quantum oscillators are around 1.4 - 1.6 [GHz]. The 
resonators here need to have high quality factors. 

Conclusion and Outlook. We have shown that it 
is appealing to simulate the dynamics of open quantum 
systems with complex environments and structured spec¬ 
tral densities (such as, the chlorosome or the examples 
given in Figure 1) by using a chain of few tens of co¬ 
herent qubits. In our previous work [23], we presented 
a detailed study on simulating the dynamics of Fenna- 
Matthews-Olson photosynthetic complex as an exam¬ 
ple of complex open quantum systems. The main fo¬ 
cus of current manuscript has been to address the limit 
that analog quantum simulators based on superconduct¬ 
ing circuits with precisely engineered quantum environ¬ 
ment may outperform exact classical computational ap¬ 
proaches, such as the HEOM approach, allowing us to 
study non-Markovian effects. Furthermore, here we have 
discussed the simulation of standard, as well as, general¬ 
ized Holstein model at finite temperature which has many 
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applications in molecular aggregates, polymers, and su¬ 
perconductivity. Using the linear-algebraic bath trans¬ 
formation, we will be able to simulate dynamics of com¬ 
plex open quantum systems with thousands of phonon 
modes. Such a simulation is as exact as numerical ap¬ 
proaches such as HEOM and definitely out of reach of 
any currently available computational device. 
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